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Abstract 



The Fermi-Pasta-Ulam a-model of harmonic oscillators with cubic anhar- 



o 

O 

■ monic interactions is studied from a statistical mechanical point of view. Sys- 
^ ■ tems of = 32 to 128 oscillators appear to be large enough to suggest 
O ■ statistical mechanical behavior. A key element has been a comparison of the 

V~j ■ maximum Lyapounov coefficient Xmax of the FPU a-model and that of the 

■ Toda lattice. For generic initial conditions, Xmax{t) is indistinguishable for 
the two models up to times that increase with decreasing energy (at fixed 
N). Then suddenly a bifurcation appears, which can be discussed in relation 
to the breakup of regular, soliton-like structures. After this bifurcation, the 
Xmax of the FPU model appears to approach a constant, while the Xmax of 
the Toda lattice appears to approach zero, consistent with its integrability. 
This suggests that for generic initial conditions the FPU a-model is chaotic 
and will therefore approach equilibrium and equipartition of energy. There 
is, however, a threshold energy density ec{N) ~ l/iV^, below which trapping 
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occurs; here the dynamics appears to be regular, sohton-Uke and the approach 
to equihbrium - if any - takes longer than observable on any available com- 
puter. Above this threshold the system appears to behave in accordance with 
statistical mechanics, exhibiting an approach to equilibrium in physically rea- 
sonable times. The initial conditions chosen by Fermi, Pasta and Ulam were 
not generic and below threshold and would have required possibly an infinite 
time to reach equilibrium. The picture obtained on the basis of Xmax suggests 
that neither the KAM nor the Nekhoroshev theorems in their present form are 
directly relevant for a discussion of the phenomenology of the FPU a-model 
presented here. 

PACS numbers(s): 05.45.+b; 05.20.-y 
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Typeset using REVT^ 



I. INTRODUCTION 



Few problems have been studied so extensively over the last decades as the one devised 
originally by E. Fermi, J. Pasta and S. Ulam (FPU) in 1954 [jl|. Their purpose was to check 
numerically that a generic but very simple non-linear many particle dynamical system would 
indeed behave for large times as a statistical mechanical system, that is it would approach 
equilibrium. In particular their purpose was to obtain the usual equipartition of energy over 
all the degrees of freedom of a system, for generic initial conditions. To their surprise, for the 
system FPU considered - a one dimensional anharmonic chain of 32 or 64 particles with fixed 
ends and in addition to harmonic, cubic (a-model) or quartic (/?-model) anharmonic forces 
between nearest neighbors - this was not observed. A variety of manifestly non-equilibrium 
and non-equipartition behaviors was seen, including quasiperiodic recurrences to the initial 
state. In fact, a behavior reminiscent of that of a dynamical system with few degrees of 
freedom was found, rather than the expected statistical mechanical behavior. The duration 
of their calculations varied between 10000 and 82500 computation steps. These results raised 
the fundamental question about the validity or at least the generally assumed applicability of 
statistical mechanics to non-linear systems of which the system considered by FPU seemed 
to be a typical example. Most of the attempts to clarify this difficulty have approached the 
problem as one in dynamical systems. These analyses have revealed many very interesting 
properties of the FPU system and uncovered a number of possible explanations for the 
resolution of the observed conflict with statistical mechanics. The classical explanations 
are: i) the survival of invariant tori in the phase space of a quasi- integrable system (KAM 
theory) ii) the existence of Zabusky and Kruskal's solitons in the KdV continuum limit 
0,^, in) the existence of an order-to-chaos transition ||^. However, we do not believe that 
the problem has as yet been entirely resolved. In particular, it is the purpose of this paper 
to try to clarify the problem from a statistical mechanical point of view. That is, we will 
try to exhibit the reasons why this apparently bona fide statistical mechanical system did 
not behave as such and, in particular, what in our opinion the significance of this apparent 
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failure is for the general validity of statistical mechanics. 

There are a number of obvious questions related to the unstatistical mechanical behavior 
observed by FPU, which all address the generic nature of their results: 

a) Was their time of integration long enough? 

b) Was their dynamical system of = 32 or 64 particles in one dimension large enough, 

i.e. possessing a sufficient number of degrees of freedom, to qualify as a statistical 
mechanical system? 

c) Were the recurrence phenomena (to within 3%) observed by FPU, transient or generic, 

i.e. possibly related to a Poincare recurrence time? 

The search for answers to these questions made the work of FPU very seminal, spawning 
many new developments and connections in the theory of nonlinear dynamical systems, 
such as the connection with continuum models based on the Korteweg - de Vries equation, 
leading to solitons heavy breathers etc., or with few degrees of freedom models like the 
Henon-Heiles and the Toda lattice |0. For a recent review we refer to 

The approach to equilibrium of the FPU - (3 model was studied extensively for various 
classes of initial conditions by Kantz et al. p| and recently by De Luca et a/. 0]. A very 
detailed picture has emerged from their work, as to the behavior of the FPU-/9 model 
in its dependence on non-equilibrium initial conditions as well as the role played by low 



frequency and high frequency mode-mode couplings during its time evolution. Two 
threshold energies were identified, but the connection of all the very interesting and detailed 
information obtained with the generic statistical mechanical behavior of the FPU-/5 model 
remains unclear so far. 

Thus, although the effort to resolve the so-called FPU problem has led to enormous 
advances in our understanding of non linear dynamical systems, it has not yet, in our 
opinion, led to a full evaluation of the statistical mechanical relevance of the FPU paradox, 
i.e. FPU's original question has not yet really been answered. 



It is commonly asserted that the KAM theorem provides the essential answer to 
FPU's observations, i.e., for sufficiently small nonlinearities and a class of initial conditions 
living on non-resonant tori, the FPU system behaves like an integrable system and is rep- 
resented by deformed tori in phase space. With increasing strength of the non-linearities, a 
progressive chaotic behavior appears, which would ultimately lead to the expected approach 



to equilibrium and equipartition |Tl 



Even though there are regular regions in phase space, the existing estimates based on the 
KAM theorem are qualitatively different from what we found, indicating that the physics 
of the FPU model is quite different from what is contained in these estimates. This makes 
us believe, on the basis of our numerical simulations, that it are the very special initial 
conditions chosen by Fermi and collaborators that make their system belong to a regular 
region of phase space. If they had chosen an initial excitation ten times larger they could 
have observed equipartition. This takes place via a disappearence with increasing initial 
amplitude of a threshold, whose iV-dependence is entirely different from that estimated in 
the KAM theory framework. Hence the source of FPU's failure is the fortuitous choice of 
initial conditions in a regular region of phase space, below this critical energy. 

Thanks to the power of modern computers, we have considerably extended the calcula- 
tions performed in the past by various authors. We have been able to reconcile different, and 
sometimes contradictory, aspects of the FPU dynamics, finding that regular regions and a 
large "chaotic sea" can coexist in phase space. The lack of equipartition in the original FPU 
experiment is not representative of a global property of phase space: apparently regular, 
soliton-like structures, similar to those of Zabusky and Kruskal, have a very long, possibly 
infinite, life-time below a stochasticity threshold, whereas, above the same threshold, they 
have only a finite life-time. 

By choosing more physically generic initial conditions, i.e random positions and mo- 
menta, a threshold energy for the onset of chaos is detected, showing a rather strong tendency 
to vanish at increasing number of degrees of freedom. Thus we have found strong evidence 
in support of the point of view that the so-called "FPU-problem" does not invalidate the 



(generic) approach to equilibrium and the vahdity of equihbrium statistical mechanics. On 
the other hand the existence of long living initial states far from equilibrium, may well have 
interesting, non trivial physical implications. 



II. MODEL AND RESULTS 



We have considered a one dimensional lattice of unit mass particles interacting via 
nearest-neighbor forces with unit harmonic coupling constant, with fixed end-points {qi = 
Qn+i = 0) and described by the Hamiltonian 



N 



H{P, q) = E 



k=l 



\pi + \{(ik+i - qkf + |(gfe+i - qtf 



(1) 



which is known as the FPU a-model. We can think of (|I|) as the first anharmonic approxi- 
mation to physical interatomic potentials of the Morse or Lennard- Jones type. 

The cubic term is obviously responsible for the energy exchange among the normal 
modes of the harmonic part of ([^). The normal mode coordinates, obtained by a standard 
orthogonal transformation, read 



Qk 



— y Qi sm 



(2) 



and diagonalize the harmonic part of (|T|). 

The Hamiltonian in these new coordinates becomes 



N 



k=l 



N 



^ ^ k'k"=i 



(3) 



where ujk = 2sin(^|^^ is the frequency of the k-th normal mode and Pk = Qk are the 
conjugated momenta. The natural unit of time with the choice of the units of Eq. (0), is given 
by the inverse of the fastest frequency of the harmonic part of (0): T^m = '^'^/^max = ^r. In 
what follows t = 1 corresponds to l/vr of this fastest linear period. 

We have numerically integrated this model by means of a very efficient third order 



bilateral symplectic algorithm |jT2| that ensures a faithful representation of a Hamiltonian 



flow and, with the adopted values of the rather large time step ranging from 0.01 to 0.1, keeps 
the total energy E constant within an average fluctuation level of AE/E ^ 10"^ without 
drift. Such a high precision in numerical integration makes the outcome of the very long runs 
that are reported in the following reliable. The coupling constant was a = 0.25. In order to 
give an idea of the computational effort that was necessary to obtain the results reported in 
this paper, let us mention that the computation of the largest Lyapunov exponents Ai at low 
energy typically required integration times in the interval 10^ - 10^ natural units of time: 
the longest runs lasted about 60 hours of CPU time on an ifP9000/735 computer (about 
the same CPU time would be necessary on a CRAY Y-MP computer). The CPU time 
amounted to about 4000 hours on the following computers: ifP9000/735, ifP9000/715, 
Sun SPARCIO, Sun SPARC5, Digital AlphaServer 2000 4/200. 

For what concerns the initial conditions, we begin by choosing single mode excitations 
as Fermi and collaborators did in their original experiment, i.e. the initial displacements of 
the particles from their equilibrium positions are given by the fundamental mode: 



and the initial momenta Pi{0) = for i = 1, . . . , N. Fermi et al. used = 32, A = 1, and 
mode number n = 1. 

Now the main question is: if we repeat the original experiment, do we have any chance 
to find equipartition or any clear tendency to it by using present day fast computers? 

To answer this question we cannot blindly make the longest integration of the trajec- 
tories of the Hamiltonian (|l]), that we can afford on a very fast computer: the absence of 
equipartition, even after a very long integration time, would not be by itself conclusive for 
the non existence of equipartition for this model. Rather, one should make an estimate 
of the equipartition time aX A = 1 and = 32, and determine the A-dependence of this 
time for larger values of A: a large initial excitation amplitude makes the anharmonicity of 
the system larger, hence the mode-mode couplings stronger, so that a faster relaxation to 
equipartition can be expected. 




(4) 
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A. Detecting energy equipartition 



A possible method to numerically detect equipartition of energy makes use of the spectral 
entropy [|13| defined by 



N/2 

S{t) = -Y.Wk{t)\nwk{t) (5) 

k=l 

where the weights Wk are given by the fraction of the total harmonic energy = ki^k + 



UJ 



"^QD in the k-th normal mode JT^ 



Ek(t) 

so that S{t) reaches its maximum value when equipartition is attained. This entropy can 
be normalized as follows 

Vit) = , (7) 

hence rj = 1 detects a "freezing" of the initial condition and = detects equipartition. 

By following the time relaxation of ri{t) we have obtained some estimates of the equipar- 
tition time at values of the initial excitation amplitude A ranging from 3 to 11. Note 
that A cannot be too large, since then the cubic part of the potential becomes repulsive and 
the phase space trajectories are "runaway" trajectories, nor can A be too small, since then 
the relaxation times become too large to be determined by our computations. However, two 
major difficulties arise with the ?]- method: i) contrary to previously investigated systems 



15| , p!6|| , the relaxation pattern of rfit) does not show clearly when equipartition sets in, so 
that the equipartition time is a fuzzy quantity, ii) if we somehow define the equipartition 
time by, for example, measuring the time needed for ri{t) to drop below a threshold value, 
say 7] = 0.1, we find that T^(e) ~ e^'^, (e = E/N is the energy density). For the FPU case 
[A = 1, e = 0.00241), we will show (see Section [II B| ) that the equipartition time is far too 



long a time for numerical tests, since it would amount to about a year's CPU time. 
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B. Phase space trapping 



We have been able to overcome these difficulties by focussing on the development of 
chaoticity in the time evolution of the system rather than on the attainment of equipartition. 

The natural way of characterizing chaotic motions is to compute the largest Lyapunov 
exponent Ai. Let us briefly remember its definition. Let 

x' = X\xK..x^) , i = l,...,N (8) 

be a given dynamical system, and denote by 

N 

e = E'^km]e, t = i,...,N (9) 

k=i 

its tangent dynamics equation, [J'^] = [dX'^/dx'^] is the Jacobian matrix of [X*], then the 
largest Lyapunov exponent Ai is defined by 



Ai = lim - In H^) . (10) 

By setting K[x{t),i{t)] ^ e JW)]i/ = ei/Ci = liHei), Ai can be formally 
expressed time average 

Ai = ;nn l|JrfrA[x(r),e(r)] (11) 



which, in practice, is evaluated by computing fl? 



where t„ = riAt (At is some time interval), up to a final time tj^, n = A/", such that Ai(t^) 
has converged to a reasonably asymptotic value [|1^]. A positive asymptotic value of Ai 
obviously detects a chaotic dynamics, whereas \i(t) ~ t^^ corresponds to a non-chaotic 
dynamics. 

We have exploited the existence of an integrable model, the Toda lattice, close to the 
FPU-a model. The Toda lattice Hamiltonian reads, in the same units as used in Eq.(|T]): 



^ rl 1 1 

^(P, q) = J2 2^fc + ^ {exp[-2a(gfc+i - g^)] + 2a{qk+i - g^) - 1} ; (13) 

and the power series expansion of its potential coincides, up to third order, with the potential 
part of (1). 

As the Toda lattice Hamiltonian describes an integrable system, the numerical compu- 
tation of the time behavior of the largest Lyapunov exponent Xi{t) must reveal the non- 
chaotic property of its dynamics. In particular, one might wonder whether some informa- 
tion can be obtained by comparing X[ [t) with a; (t) for the same initial conditions, i.e. 
{gi(0),pi(0)}, i = 1, . . . , N. This turns indeed out to be the case. In Fig. 1 we report, as an 
example, A^^'^(t) and Xl°''"'{t) in the case = 32 and e = 0.0217 for an initial condition with 
y4 = 3 in Eq. @). Up until t ^ 2 ■ 10^ these two functions are so close to each other that the 
two are virtually indistinguishable. Then - suddenly - they separate at t > 4 ■ 10^: A^°''''(t) 
continues its decay toward zero, whereas Xi''" (t) tends to converge to a non- vanishing value. 
This makes possible to define clearly what a trapping time in a regular region of phase space 
is; moreover, its numerical determination is unambiguous, as it can be deduced by simply 
looking at Fig. 1. We attribute this dramatic difference to the un-trapping of the FPU 
system from its regular region in phase space by escaping to the chaotic component of its 
phase space. 

The peculiar behavior of Xi^" (t) suggests therefore that non-integrable motions, origi- 
nated by one-mode initial excitations, after a transient, possibly long, trapping in a regular 
region of phase space, enter its chaotic component. The chaotic component of phase space, 
by the Poincare- Fermi theorem is connected, so that we may well expect that equipar- 
tition is eventually attained on a finite, albeit possibly very long, time scale. 

The precise nature of the trapping mechanism is unclear to us. One could either think 
that the trajectories stick close to a KAM torus during the trapping time, or that - for 
example - they are geodesies of a bumpy manifold to which an A^-torus has been differentiably 
glued by a tiny "bridge"; in that case a geodesic - originating at any point of the A^-torus - 
is locally stable until it finds a path to escape from the torus. Another possible mechanism 
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of trapping and escape is discussed in the Discussion Section under point 5, in the hght of 
the data discussed in the sequel of the present Section. 

In conclusion, once we observe that a trajectory enters the chaotic component of phase 
space, we are allowed to think that equipartition will be eventually attained. In Fig. 2 the 
relaxation times to equipartition TE{e) and the trapping times Tj'(e), evaluated for the FPU- 
a model, are shown. The results refer to = 32 and initial conditions given by Eq.(^). The 
initial excitation amplitudes A range from 3 to 11. The equipartition times appear to be 
between one and two orders of magnitude larger than the trapping times. If we extrapolate 
pO| the equipartition time to the case of FPU, we find ~ 4 ■ 10^'^ (as is indicated in Fig. 
2 by an asterisk). 

Since we can consider the FPU a-model as a third order expansion of a Lennard- Jones 
potential around its minimum, we can tentatively extrapolate the equipartition times re- 
ported in Fig. 2 to a macroscopic system in three dimensions. As an example, we can 
roughly estimate what the physical equipartition time could be for a classical Xenon crystal 
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at zero temperature and only one normal mode initially excited. At 6 — 0.01, ~ 10 
proper times (see Fig. 2) corresponds to about 10~^s, using as proper time, obtained with 
standard Lennard- Jones parameters for Xenon, 2.4 ■ 10~^^s. 



C. Stochasticity threshold 

We have evaluated the trapping times T^{e,N) for the FPU-a model at different values 
of both the energy density e and of the number of degrees of freedom A^. When A^ was 
varied, we kept the wavelength of the initial excitation constant (i.e. n in Eq.(^ is taken 
proportional to A^: n = 1 aX N = 32, n = 2 at A^ = 64 and n = 4 at A^ = 128), and we excited 
only one mode at t = 0. The results are reported in Fig. 3. With decreasing e, first tends 
to increase monotonically, then, abruptly, it displays an apparently divergent behavior. In 
fact, when reaching critical threshold values, a very small change in e suddenly gives an 
extremely steep increase in r^. This very steep increase of with decreasing e suggests at 
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least a very narrow bottleneck in phase space, through which the system can only escape 
with great difficulty. We assume that this bottleneck is not an insurmountable barrier, and 
for that reason, as well as to conform with previous use, we will call it a threshold. This is 
consistent with the results of Fig. 4, where the largest Lyapunov exponents are reported for 
the same cases as in Fig. 3. We have used a cut-off of the integration time at t = 4.3 ■ 10^. 
After such a long time - when e was smaller than the threshold value - no separation between 
Xj^'^^it) and Af^^(t) has been observed: the values of Af^^ at t = 4.3 ■ 10^ (endpoints of 
broken lines in Fig. 4, marked by arrows) are taken as upper bounds for the FPU-Lyapunov 
exponents and t = 4.3 ■ 10*^ is taken as a lower bound for the trapping time of the FPU- 
phase point (endpoints of the broken lines in Fig. 3, marked by arrows). Both Ai(e, A^) 
(from here on, by Ai we mean Af'^'^) and T^{e,N) strongly suggest the existence of an 
A^-dependent threshold value of the energy (density), above which the motion is chaotic 
and below which the trajectories appear to be trapped in a regular region of phase space 
[Stochasticity Threshold (ST)]. 

The following approximate relationship between A and e holds: e ~ 0.00241^4^ therefore 
the threshold amplitudes Ac leading to chaos are: Ac{N = 32) ~ 1.62, Ac{N = 64) ~ 1.42, 
Ac{N = 128) ~ 1.28, respectively. 

The dotted line at e = 0.00241 corresponds to the initial excitation amplitude A = 1 of 
FPU's original paper. We see that Fermi and coworkers chose an initial condition well below 



this ST pO[. Figure 3 shows that if they had taken ten times larger amplitude, they would 
have observed equipartition during the integration time they used. This appears to us to be 
the explanation of the lack of statistical mechanical behavior observed in the original FPU 
numerical experiment. 

A few more comments on these results: i) Fig. 3 suggests that with increasing A^, 
Tj, (e, A^) probably becomes less dependent on A^ [in fact the values of (e, A^ = 64) are 
closer to those of T^{e, N = 128) than to those of T^{e, N = 32)]; ii) there is a narrow energy 
density interval where and Ai oscillate, suggesting the transition in phase space from 
chaotic behavior at large e to regular behavior at small e; in) the ST shows a weak but clear 
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dependence on (provided the wavelength of the initial excitation is kept constant). This 
is not surprising because if we combine two identical systems - each with the same initial 
excitation - the composite system will have new low frequency modes that are absent in the 
separate subsystems, which can be expected to facilitate the mechanism of energy exchange 
among the normal modes. 

As far as we are aware, the direct evidence given in Fig. 4 of the existence of a ST at 
N ^ 2, obtained through the behavior of the largest Lyapounov exponent, was never found 
before in nonlinear Hamiltonian systems. An indirect suggestion of its existence (through a 
"freezing" of the decay of the spectral entropy) has been given in for the FPU /3-model. 

D. Coexistence of order and chaos 

A question now arises: does such a threshold refer to a global property of the constant 
energy surface or is it, rather, a local property of sensitive to the initial condition? 

In order to answer this question we have considered also the following initial conditions 
at A^ = 32: i) A ^ [the fundamental mode - Eq. (^) - is excited with amplitudes 
ranging from A = 0.7 to A = 5.5 ] and Pi{0) ^ {i = 1,...,N) are gaussian random 
numbers with zero mean and standard deviation ("temperature") equal to 0.001; ii) A = 0, 
qi{0) = {i = 1,...,N) and Pi{0) {i = 1,...,N) are randomly chosen according to a 
gaussian distribution. 

In the first case, a fraction of the energy in the initial excitation is given to all the normal 
modes of the system: in so doing, we are displacing the starting point on S^, proportional 
to the magnitude of the standard deviation of the noisy component of the initial excitation. 
This is intended to provide some information about the extension of the regular region(s) 
in phase space. If all, or almost all, the energy of the initial excitation is concentrated 
in one or a few modes, then we are dealing with a non-equilibrium initial condition; the 
existence of a ST entails then that a system, prepared in a non-equilibrium state below such 
a threshold, will appear never to attain equipartition of the energy of the initial excitation. 
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The noisy component has a self evident physical meaning related to the impossibility of 
preparing any physical system in a perfectly ordered initial state: at non-zero temperature 
some randomness in the initial conditions is unavoidable. 

The second choice - completely random initial conditions - mimicks a physical situation 
that corresponds to a crystal prepared at an assigned temperature (i.e. mean kinetic energy 
per degree of freedom) at thermal equilibrium. 

In Fig. 5 the effect of changing the initial conditions according to the above prescriptions 
is shown. Here = 32, the squares refer to A ^ and no random component, the asterisks 
refer to only random initial excitations and the star-like squares refer to A 7^ plus a random 
component. 

In each case a threshold energy (or equivalently energy density, since N is fixed) is found. 
At e > 0.01, the uncertainties in the determination of Ai are of the order of the size of the 
symbols used; at e < 0.01, an estimate is difficult because of unpredictable fluctuations of 
Ai(t) that could be reduced only by prohibitively long integration times. Nevertheless, the 
information given by Fig. 5 is unambiguous. Down to e ~ 10~^ all the values of Ai(e), 
obtained with different initial conditions, crowd along the same line; below e ~ 10^^ the 
values of Ai, obtained with different initial conditions, separate and exibit a coexistence of 
regular and chaotic regions on the constant energy surface in phase space, whose details 
depend on the initial conditions. Below e ~ 1.8 ■ 10~^, even with only random initial 
conditions, the motions are regular (star-like squares). 

This tells us that the phase space undergoes some important structural change as a 
function of the energy, in analogy with what is observed in two-degrees-of-freedom systems. 



like the Henon-Heiles model where fully developed chaos, a coexistence of regular and 
chaotic regions of phase space, or only regular trajectories are observed, depending on the 
value of the energy. 
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E. N-dependence of the stochasticity threshold 



An important question is the stability or instability of the stochasticity threshold with 
A^. To this end we have numerically determined Ai(e, A^) at A^ = 8, 16, 32, 64 always choosing 
random initial conditions, i.e. {qi{0) = 0,pj(0) = rj}, i = 1, . . . , N, with rj gauss-distributed 
random numbers with zero mean and variance \/2e (after the assignement of the random 
values Tj, the momenta are adjusted, by rescaling them, in order to obtain a fixed initial e). 

In Fig. 6 the outcome of these computations is reported. The endpoints of the broken 
lines have the same meaning as already discussed above. At large e there is a tendency of 
all the sets of points to join. This fact is most evident for A^ = 32 and A^ = 64 which have 
a line segment in common and then separate at small e: the larger A^, the smaller the e at 
which the separation occurs. 

At each A^, we take as rough estimate of the stochasticity threshold ec, the value of e at 
the midpoint between the two lowest points on each curve, because the lowest point (marked 
by an arrow) is presumably below threshold. 

Figure 7 then shows that for these threshold values Cc plotted vs. N holds: edN) ~ l/N"^. 
This result is interesting for the following reasons: 

i) the threshold values vanish sufficiently fast with increasing A^, that the existence of 

regular regions of phase space below Sc does not constitute a problem for equilibrium 
statistical mechanics. In fact, it appears that A^ = 32, . . . , 128 is not too small to 
obtain indications, if not "confirmation" , of the statistical mechanical behavior of the 
FPU-a system. 

ii) both the values of ec and the A^-dependence €c{N) ~ 1/A^^ can hardly be explained on the 

basis of the best available estimate of the perturbation amplitudes for which a positive 
measure of the regular KAM regions in phase space exists: /i < /ic ~ aexp(— 6A^ln A^) 
^], where /i measures the relative strength of the anharmonic to the harmonic part of 



a given Hamiltonian (/i depends on e), /ic is a threshold value, a and b are constants. 
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There are also power-law estimates for fic{N), i.e. /ic(A^) ^ N ^, obtained in the 
context of KAM theory [0], however - at present - 6 is still very large: 6 ~ 160. 

iii) the bounds on the threshold value < Cc < 2c?lv fo^^d by Enz et al. for the 

FPU-a model predict an A^-dependence which is in much better agreement with 
our results than the exponential drop with mentioned in the previous point. 



III. DISCUSSION 

We conclude with the following remarks. 

1. It is worth mentioning that recently the existence of an equipartition threshold vanish- 
ing at increasing A^ has also been reported for the FPU-/? model PJ2^J2^ . It is inappropriate 



to compare these results quantitatively with ours: our model is different and the evidence 
for a stochasticity threshold for the FPU /5-model has been obtained indirectly, through the 
opening of a local trap in phase space, which prevented equipartition. However, we can 
say that both results are in qualitative agreement: threshold effects, either concerning tran- 
sient dynamics to equipartition of non-equilibrium initial conditions (studied through spec- 
tral entropy) or concerning the dynamics of equilibrium initial conditions [studied through 
Ai(e, A^)], vanish at increasing A^. 

A qualitative agreement about the vanishing with A^ of the critical energy to get chaos 
is reported in a recent paper on the FPU a-model p7| . 

The question of how to explain the existence and the 1/A^^ dependence of the stochasticity 
threshold reported here, remains open. 

2. In this paper we report the existence of two interesting phenomena among others: 
the apparent existence of regular regions in the phase space of a non-integrable Hamiltonian 
system as is the FPU-a model and the existence of almost regular regions of phase space 
where the trajectories are trapped during long but finite times. These phenomena are 
reminiscent of the KAM and Nekhoroshev theorems respectively. We have already 
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discussed throughout the paper why our results disagree with the present-day quantitative 
predictions of the KAM theorem. Similarly, the Nekhoroshev theorem does not seem to be 
able to provide an explanation of the observed phenomenology as well. In fact, without 
entering into the details of how a thorough comparison with our results could be made, 
Nekhoroshev's estimate of the lower bound of the "trapping" time T of a trajectory close 
to its initial condition is T ^ exp(c//i)'^*-^'', where c is a constant, /i has the same meaning 
as above, and 7(A^) ~ 1/8A^ is the optimal A^-dependence of the exponent. Hence 



for = 32,64, 128, as is the case of the results of Fig. 3 for T^{e,N), the "Nekhoroshev 
trapping" time T is 0{1) independently of fi (thus of e). Therefore, it appears that the 
physical mechanisms responsible for finite time trapping in phase space, that are behind the 
Nekhoroshev theorem and our numeric phenomenology, respectively, are different. Thus we 
emphasize the difference between, on the one hand, the approach to infinite time trapping 
(KAM theorem) and to finite time trapping (Nekhoroshev theorem) based on the description 
of the persistence of certain local properties of regular regions of phase space, and, on the 
other hand, our chaos- related approach, based on the comparison of Xi{t) for the FPU and 
Toda lattices (Fig. 1). The behavior of Xi{t) suggests that the sudden escape from the 
regular region occurs as if the trajectory would eventually find a "hole" in its boundary. 

3. It is not out of place to note that while the ST drops to zero as A^ — > oo in the FPU 
a-model, there is another transition phenomenon that does not suffer this A^-dependence: 
it is the Strong Stochasticity Threshold (SST) that concerns a transition from weak to 
strong chaos p!5| , p!6| which has been discovered in the FPU /3-model and in the lattice ip^ 
model. Unfortunately, this SST cannot be investigated in the FPU-a model because the 
cubic potential prevents working at too large energy. We note, however, that the stability 
with A^ of the SST makes the SST not only a dynamical phenomenon related to dynamical 
chaos, but also of potential interest to statistical mechanical phenomena. 

4. The existence of special initial conditions that create non-trivial dynamical behavior, 
though irrelevant for equilibrium statistical mechanics because of their negligible measure, 
could be physically relevant for transient non- equilibrium statistical mechanics if we can 
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conceive an operational method to prepare a real system in such a special initial state (see 
for example pO| ). 

5. Let us comment about the classical explanation of the non-statistical mechanical 
behavior of the FPU model proposed by Zabusky et al. based on the existence of 

soliton solutions of the Korteweg-de Vries (KdV) equation, derived as a special continuum 
limit of the FPU a-model. Their experiments were carried out over much shorter time 
scales than those that are possible nowadays. As we have seen above, thanks to very long 
numerical integrations, it turns out that, below a threshold, regular regions of phase space 
can coexist with chaotic ones. In the light of our results, we can thus assert that the KdV 
soliton solutions belong to the regular region of phase space of the FPU system from which 
- after sufficiently long time - escape will occur to the chaotic component of phase space. 

In Fig. 8 the patterns of the displacements as a function of position i are shown at 
different times in the case = 32, A = 3. For times t below the trapping time tt — 4 ■ 10^, 
they are looking as regular structures, apparently composed of a superposition of a small 
number of waves, which display a trapping and a soliton-like recurrence (see Fig. 8b), sim- 
ilar to that observed by Zabusky and Kruskal |^ and Tuck and Menzel [0, although we 
have fixed boundary conditions, as in FPU's original paper, rather than periodic boundary 
conditions as considered by Zabusky and Kruskal. For t > tt we observed a gradual un- 
trapping or decay of these regular features due to more and more complicated structures 
consisting of a superposition of an increasing number of different waves (radiated by the 
decaying solitons ^). Finally, for t > 10^, an almost noisy pattern is attained (see Fig. 8c) 
as well as energy equipartition (detected through the spectral entropy). This demonstrates 
clearly the compatibility of the existence of trapping, i.e. very long lived regular solutions, 
and the attainment of equipartition, albeit after possibly very long times. It is still an open 
question whether, at fixed N, below the Stochasticity Threshold, i.e. at energy density such 
that the largest Lyapunov exponent seems to vanish for any initial condition, the lifetime of 
regular solutions might actually diverge. 

However, we have found that this threshold energy density drops to zero as ~ l/N"^. 
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Therefore it appears to us that at sufficiently large N the existence of KdV sohtons does 
not hinder a good statistical mechanical behavior of the FPU system. 
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FIGURES 

FIG. 1. X^''" (t) (solid triangles) and A^°'*"(t) (squares) are plotted vs. time for = 32 and 
energy density e = 0.0217 (which corresponds to an initial excitation amplitude ^ = 3 in Eq. (^). 
Dividing t on the horizontal axis by vr gives the time in units of the fastest period of the harmonic 
part of the chain, Tmin = tt. 

FIG. 2. The relaxation times to equipartition (full squares) and the trapping times (open 
squares) are plotted (for the FPU-a model) vs. the energy density e for = 32 and the ini- 
tial conditions of Eq. (^). The initial excitation amplitudes considered are, from left to right, 
A = 2,3, 4, 5, 8, 10, 11, respectively. The asterisk represents the extrapolation of the equipartition 
time to the case A = 1 (FPU's original paper [pO| ). 

FIG. 3. The trapping times tt{£,N) at different values of energy density e (i.e. at different 
values of the initial excitation amplitudes A), are reported. Open squares refer to the case = 32 
(A ranges from 1.6 to 11), solid triangles refer to = 64 {A ranges from 1.4 to 10), open circles 
refer to A^ = 128 {A ranges from 1.25 to 9), respectively. The endpoints of the broken lines are 
lower bounds for the trapping time (the cut-off of the integration time is at i = 4.3 • 10^). The 
dotted vertical line at e = 0.00241 corresponds to the initial excitation amplitude ^ = 1 of the 
FPU's original paper. 

FIG. 4. The largest Lyapounov exponents Ai(e, A^) are shown for different values of the energy 
density e and a sine wave initially excited [Eq.(4)]. Open squares refer to A^ = 32 and n = 1, full 
triangles to A^ = 64 and n = 2, open circles to A^ = 128 and n = 4, respectively. The endpoints 
of the broken lines, marked by arrows, are upper bounds for the FPU- Lyapounov exponents at 
t = 4.3 • 10^ (cut-off of the integration time). The dotted vertical line at e = 0.00241 corresponds 
to ^ = 1 (FPU's original paper). 
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FIG. 5. The largest Lyapounov exponents Ai(e) are plotted, for N = 32, at different values 
of e and different initial conditions. Squares refer to the cases A ^ and no random initial 
conditions, asterisks refer to the cases A = and random initial conditions, star-like squares to 
A^ and random initial conditions with zero mean and standard deviation 0.001. The endpoints 
of the broken lines, marked by arrows, are upper bounds for the FPU- Lyapounov exponents at 
t = 4.3 • 10^ (cut-off of the integration time). 

FIG. 6. The largest Lyapounov exponents Ai(e, A^) are plotted vs. the energy density e, for 
different values of N. Random initial conditions are chosen. Star-like polygons refer to iV = 8, 
crosses to N = 16, asterisks to N = 32, star-like squares to AT = 64, respectively. The endpoints 
of the broken lines, marked by arrows, are upper bounds for the FPU- Lyapounov exponents at 
i = 4.3 • 10^ (cut-off of the integration time). 

FIG. 7. The values of the stochasticity thresholds Cc are plotted for different values of N (for 
the estimate of Sc, see text). 

FIG. 8. Snapshots of the displacements Qi of 32 FPU-a particles as a function of position 
i along the chain shown at different times. Here the initial amplitude is ^4 = 3 - as in Fig. 1 - 
for a sine wave [see Eq.(4)]. (a) ordered configurations at t = (I), 10^ (II) and 5 • 10^ (HI); (b) 
ordered configurations at t = 10^ (I), 2 • 10^ (II) and 4.2 • 10^ (HI), note the almost recurrence at 
t = 2- 10^; (c) t = 1.6 ■ 10^ (I), 5.4 • 10^ (II) and 2 • 10*^ (III). Note the increasing disorder of the 
configurations in (c), corresponding to the onset of a chaotic dynamics in the system, due to the 
breakup of regular soliton-like structures, and to its approach to equipartition. 
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